Semiclassical quantization of an A^-particle Bose-Hubbard model 
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A semiclassical Bohr-Sommerfeld approximation is derived for an A'^-particle, two-mode Bose- 
Hubbard system modeling a Bose-Einstein condensate in a double-well potential. This semiclassical 
description is based on the 'classical' dynamics of the mean-field Gross-Pitaevskii equation and is 
expected to be valid for large A'^. We demonstrate the possibility to reconstruct quantum properties 
of the A-particle system from the mean-field dynamics. The resulting semiclassical eigenvalues and 
eigenstates are found to be in very good agreement with the exact ones, even for small values of N , 
both for subcritical and supercritical particle interaction strength where tunneling has to be taken 
into account. 

PACS numbers: 03.65.W, 03.65.Sq, 03.75.Lm 



I. INTRODUCTION 

Even for weakly interacting particles, a full many-body 
treatment of Bose-Einstein condensates (BEG) is only 
possible for a small number TV of particles. Most often 
a mean-field approximation is used, which describes the 
system quite well for large N at low temperatures. In this 
mean-field approach, the bosonic field operators are re- 
placed by c-numbers, the condensate wavefunctions. This 
constitutes a classicalization and therefore the result of 
the mean-field approximation, the Gross-Pitaevskii equa- 
tion (GPE), is often denoted as 'classical', despite of the 
fact that the GPE is manifestly quantum, i.e. it reduces 
to the usual linear Schrodinger equation for vanishing in- 
terparticle interaction. Therefore, in order to avoid mis- 
understanding, this inversion of the second quantization 
may be called a second classicalization. 

In a number of recent papers, consequences of the 
classical nature of the mean-field approximation are dis- 
cussed and semiclassical aspects are introduced. For a 
two-mode Bose-Hubbard model, Anglin and Vardi [H, Q 
consider equations of motion which go beyond the stan- 
dard mean-field theory by including higher terms in the 
Heisenberg equations of motion. The classical-quantum 
correspondence has been studied in terms of phase space 
(Husimi) distributions Q for such systems. Mossmann 
and Jung [3] demonstrate for a three-mode Bose-Hubbard 
model that the organization of the A^-particle eigenstates 
follows the underlying classical, i.e. mean-field, dynamics. 
A generalized Landau-Zener formula for the mean-field 
description of interacting BEGs in a two-mode system 
has been derived by studying the many particle system 
0. In the commutability between the classical and 
the adiabatic limit for the same system is studied and 
first steps towards a semiclassical treatment of the prob- 
lem are reported. 

The purpose of the present paper is to show that the 
mean-field model is not only capable to approximate the 



interacting A^-particle system in the limit of large N 
and to allow for an interpretation of the organization 
of the iV-particle eigenvalues and eigenstates, but can 
also be used to reconstruct approximately the individ- 
ual eigenvalues and eigenstates in a semiclassical Bohr- 
Sommerfeld (or EBK) manner with astounding accuracy 
even for a small number of particles. This will be demon- 
strated for TV bosonic particles in a two-mode system, a 
many-particle Bose-Hubbard Hamiltonian, describing for 
example the low-energy dynamics of a BEG in a (possi- 
bly asymmetric) double-well potential. This system is 
related to a classical non-rigid pendulum in the mean- 
field approximation (see, e.g., 0] and references therein). 

Both the many particle model and its classical ver- 
sion - which is often denoted as the discrete self-trapping 
equation for reasons which will become obvious in the 
following - have been studied for decades in different 
context (see 0)- A detailed semiclassical analysis is, 
however, missing up to now. Previous work in this di- 
rection concentrated on the symmetric case, where the 
permutational symmetry of the system with respect to 
an interchange of the two modes simplifies an analysis. 
Semiclassical expressions for the tunneling splittings of 
the eigenvalues have been derived [1, in context of 

Tl| for 
3l for a 



the spin-system in eqn. ^ below (see also [l^ 



a perturbative treatment of the splittings and 
detailed analysis of the quantum spectrum). 

In the following we will first give a short overview of 
the many particle description of the model, the celebrated 
mean-field approximation and their correspondence. Af- 
terwards we focus on the question to which extent many 
particle properties can be extracted from the mean-field 
system by an inversion of this 'classical' approximation in 
a semiclassical way using the EBK-quantization method 



II. TWO-MODE BOSE-HUBBARD MODEL 
AND MEAN-FIELD APPROXIMATION 
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In a two-mode approximation at low temperatures a 
BEG in a double-well potential can be described by a 
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FIG. 1: (Color online) Phase space portrait of the mean-field 
Hamiltonian Ti.{p, q) in ([5| for v — 1 and e — —0.5 in the 
subcritical {g — —1/Ns, top) and supercritical {g = —3/Ns, 
bottom) region for = 10 and h = 1. 



second quantized many particle Hamiltonian of Bose- 
Hubbard type: 



H 



;{hi — fi2) + v{a\a2 + a\ai) + g {n\ + . (1) 



Here hj , are bosonic particle annihilation and creation 

operators for the jth mode and iij = cijUj are the mode 
number operators. The mode energies are ±e, v is the 
coupling constant and g is the strength of the onsite 
interaction. In order to simplify the discussion, we as- 
sume here that v is positive and g is negative [14| . The 
Hamiltonian ([T|) commutes with the total number oper- 
ator N = fii + n2 and the number N of particles, the 
eigenvalue of N, is conserved. For a given value of N, 
we have N + 1 eigenvalues of the Hamiltonian ^ . Alter- 
natively, the system can be described in the Schwinger 
representation by transformation to angular momentum 
operators Jx = {a\a2 + O'l'^i)/^ , Jy = {a\a2 — a^ai) / 2\ ^ 
Jz = {o\ai — 0)20,2) The Hamiltonian ^ then takes 
the form 



(2) 



where the total angular momentum is J = N /2. 

The celebrated mean-field description can be most 
easily formulated as a replacement of operators by c- 



numbers a. 



mute in contrast to the quantum mechanical operators. 



Since the c-numbers com- 



the transition quantum — > classic and vice versa is not 
one-to-one. To avoid ambiguities one has to replace sym- 
metrized products of the operators by the corresponding 
products of c-numbers. Therefore we will start on the 
many particle side with a symmetrized Bose-Hubbard 
Hamiltonian in the following, where the fij are replaced 
by n| = (d]dj -I- dja])/2 (see also Q). This symmetriza- 
tion affects only the nonlinear term in ([T]) and the sym- 
metrized H is related to ([1]) by an additive constant term 
depending only on N. Note that thus the number opera- 
tor iV = ni+n2 = nf-l-nl— i is replaced by |?/'iP + |V'2p — 1 
and therefore the mean-field wavefunction is normalized 
as |V'i|2 + |V2p=^ + l=A^.. 

The mean-field time evolution is given by the two level 
nonlinear Schrodinger equation, resp. GPE, 



^1 



, (3) 



where "01 and ip2 are the amplitudes of the two conden- 
sate modes. 

The Schrodinger equation, linear or nonlinear, has the 
canonical structure of classical dynamics [l^, [l^, [l3| : The 
time dependence of the complex valued mean-field am- 
plitudes can be written as canonical equations of motion 



i^-r^ = TTTT and = (4) 



dt 



an 



with a Hamiltonian function Ti 

|4 



u (■(/']' ■02 + ■'/'iV'2) + ffdV'il'* + |V'2|'*)- The conservation of 
the particle number allows a reduction of the dynamics to 
an effectively one-dimensional Hamiltonian evolution by 
an amplitude-phase decomposition ipj = y^nj + 1/2 e'^^' 
in terms of the canonical coordinate q = {qi — 'Z2)/2, an 
angle, and the (angular)momentum p = (ni — 712)^., with 
the Hamiltonian function 



nip,q) 



p 
h 



cos(2(7) - 



2^^ h?' 



(5) 



where Ng is the normalization of "tp. Introducing the new 
variables the canonical equations of motion ^ obtain 
their usual appearance q = dH/dp and p = —dH/dq: 



p=2i,\/iV|-^sin(2<7) 



e 



cos(2(7) -I- g 



(6) 
(7) 



This describes the classical dynamics of a non-rigid pen- 
dulum where the phase space is finite, —Ngh < p < Nsh, 
< 9 < TT, if the lines q ~ and q — tt are identified. 

One of the prominent features of the two-mode sys- 
tem is the self-trapping effect, which leads to the emer- 
gence of additional stationary states favoring one of the 
wells above a critical particle interaction strength. For 
a discussion of the relation between mean-field and iV- 
particlc behavior see, e.g., [H, [ll] and references therein 
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and [1^ for its control by external driving fields. The 
self-trapping transition occurs at g = —v/Ng and is con- 
nected to a bifurcation of the stationary states, the fixed 
points of the Hamiltonian ([5]), in the mean-field approx- 
imation. Figure [T] shows phase space portraits of Ti.{p, q) 
for sub- and supercritical particle interaction strength. In 
the subcritical region one has a maximum, at g ~ 
and a minimum, E~ , at q ~ 7r/2. For the symmetric 
case e = 0, both are located at p = 0, which means that 
the population in both wells is the same. In the super- 
critical region the minimum bifurcates into two minima, 
E^, and a saddle point, E'^^^^^ > E^. Even for the case 
e = the two minima are located at a finite value of the 
population imbalance. In these stationary states the con- 
densate mainly populates one of the wells, which leads to 
the name self-trapping. In phase space, the regions with 
oscillations around one of the two minima are separated 
by a separatrix passing through the saddle point. The 
period of the separatrix motion is infinite. 

Figure [2] shows an example of the many particle eigen- 
values En in dependence on e for a subcritical interac- 
tion strength. The pattern of eigenvalues varies smoothly 
with e and is bounded by the stationary mean-field ener- 
gies shown as red curves. Because of the symmetry of the 
spectrum for s —>■ — e the exact spectrum is only shown 
for e < 0, whereas for e > the semiclassical eigenvalues 
are shown as discussed below. Figure 3] shows a similar 
plot in the supercritical region. Here we observe a net 
of avoided crossings clearly organized by a skeleton pro- 
vided by the stationary mean-field energies, as reported 
before by several authors d, Hi], . Again, for e > the 
semiclassical eigenvalues are shown, which closely agree 
with the quantum ones in all details. 

The mean-field eigenenergies (red curves) show a swal- 
lowtail structure which forms a caustic of the multi- 
particle eigenvalue curves in the limit N ^ oo. To il- 
lustrate this issue, one can calculate the level density 
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FIG. 2: (Color online) Many particle energies En and mean- 
field eigenenergies TCstat (- -) as a function of the onsite energy 
e in the subcritical region {g — — 0.5 /Ns) for v = 1 and — 
10 particles (ft = 1). The spectrum is organized by the mean- 
field eigenenergies Hstat (- -)• The exact quantum eigenvalues 
shown for £ < (- ) are almost exactly reproduced by the 
semiclassical ones, E^'^, shown for e > ( — ). 




E E 

FIG. 3: Level density g{E) of the many particle system in 
comparison to the mean-field energies for A'' = 1500 particles 
for v = l, g = ~3/Ns and e = 0.5, 1, 1.5 (ft = 1). 



p(E) (normalized to unity) as a function of the energy 
[la l . Figure [3] shows a histogram of the level density for 
N = 1500 particles and different values of e. The mean- 
field swallowtail curve between the cusps manifests itself 
as a peak in the density of the many particle energies. In 
the limit N oo the density g{E) approaches a smooth 
curve and the peak develops into a singularity. At the po- 
sitions of the other mean-field eigenenergies one observes 
finite steps. Indeed the quantum level densities shown in 
Fig. [3] for a large value of N are directly related to the 
classical period T of motion by T = dS/dE, where S is 
the classical action, which we will focus on in more detail 
in the following. The height of the steps in the density 
plots are simply given by the period of harmonic oscil- 
lation in the vicinity of the extrema and the singularity 
corresponds to the separatrix motion. 




FIG. 4: (Color online) Many particle energies En and mean- 
field eigenenergies Hstat (- -) as a function of the onsite en- 
ergy e in the supercritical region {g = —3/Ns) for v = 1 and 
A*' = 10 particles (ft = 1). The spectrum is organized by 
the mean-field eigenenergies Tistat (- -)• The exact quantum 
eigenvalues shown £ < ( ) are almost exactly reproduced 
by the semiclassical ones, E^ , shown for e > ( — ). 
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III. SEMICLASSICAL QUANTIZATION 
A. The classical action 

The most important ingredient of a semiclassical quan- 
tization is the action S{E), i.e. the phase space area en- 
closed by the directed curve Ti{p,q) = E. The action 
S{E) increases with E from zero at the minimum energy 
of T-L{p, q) to 2TTNsh, the total available phase space area, 
at the maximum energy of Ti.{p, q). 

For the generalized pendulum Hamiltonian (O, one 
can express the position variable q uniquely as a func- 
tion of p and E and write down the action in momentum 
space in the form S{E) ~ § q{p, E) dp. It is convenient 
psi. [2^ to introduce two functions U+{p) ~ H{p,0) and 
U-.{p) — 'H{PtTt/2), which join smoothly at p = zthNg 
and act as a potential for the variable p. The classically 
allowed energy region is given by U-{p) < E < U+{p) 
as illustrated in Fig. [5] in the sub- and supercritical re- 
gions. For a given energy E the classical turning points 
pj^ (with p_ < p+) are determined by U-{p±) = E ot 
U+{p±) ~ E. One has to distinguish three basic types of 
motion and, with S = J^* q{p, E) dp, we find: 

a) Orbits encircling a minimum of Ti.{p,q). The clas- 
sical turning points both lie on [/_ and we have 
S{E) = 7r{p+-p^)~2S. 

b) Orbits encircling a maximum of T-l{p,q). The clas- 
sical turning points both lie on C/+ and wc have 
S(E) = 2nNsh-2S. 

c) Rotor orbits extending over all angles q. We can 
find p_ on U+ and p+ on U- with S{E) = TT{Nsfi + 
p+) — 2S or p_ on [/_ and p+ on [/+ with S{E) = 
T:{Nsh-p^) - 25. 



B. Energy quantization 

In the case of a single classically accessible region, 
where there are two real turning points for any energy 
E, the semiclassical quantization condition is given by 

S{E) = h{n+\) , 71 = 0, 1, ...,7V. (8) 

This simple case is always met in the subcritical region 
\g\ < v/Ng. A numerical solution of ([5]) determines the 
semiclassical energies En, n = 0, . . . , N, where the total 
available phase space area, < S{E) < hNg, restricts the 
number of semiclassical eigenvalues to A^^ , exactly as the 
quantum ones. The resulting semiclassical eigenvalues 
shown in Fig. [2] for = 10 particles {g = —0.5/Ns, 
V = 1, h = 1) are in excellent agreement with the exact 
quantum ones. 

It should be pointed out, that in the noninteracting 
case, g = 0, the action S{E) is a linear function of the 
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FIG. 5: (Color online) The potentials U-{p) (— ) and U+{p) 
(- -), for A'' = 10 particles, v — 1 and e = 0.6 in the subcritical 
{g — ~0.6/Ns, top) and supercritical {g = — bottom) 
region, with ft = 1. 



energy E, and the semiclassical eigenvalues agree with 
the exact ones 



E„ = Ve'^ +v^{2n- N) , n = 0, 1, . . . , A^ . (9) 

This can be easily understood by recognizing that in this 
case the Hamiltonian ([1]) describes nothing but a sys- 
tem of two coupled harmonic oscillators, which can be 
transformed to two uncoupled ones by introducing nor- 
mal coordinates. It may also be of interest to note that 
(for g = and A^ = 2 or 3) the classical analog © to the 
quantum Hamiltonian ([1]) has been suggested many years 
ago by Miller and coworkers and applied in a semiclassi- 
cal description of electronic transitions in atom-molecule 
collisions [2^, [2^. 

The supercritical region \g\ > v/Ns is more compli- 
cated. Here the energy surface has two minima, hence 
the potential function U-{p) has two minima as well, 
separated by a potential barrier. In this case one has to 
distinguish different regions of the energy. For energies 
below the upper minimum (region I in Fig. [5]), the quan- 
tization can be carried out like in the subcritical case by 
equation ([8]). For energies between the upper minimum 
and the barrier E'barr (regions II in Fig. [5]), there are four 
real turning points p^^ < p'j^ < p[^^ < p!^'. In this case 
one has to consider tunneling through the barrier. The 
semiclassical quantization condition can be achieved by 
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a more elaborate expression (see also [28l [29|): 

Vl + K^ cos{Si + Sr- S4,) ^ cos(S'^ - 5,. + Se) (10) 

where Si and Sr are the action integrals in the left resp. 
right region in Fig. [5] (note that also here one has to 
distinguish the different cases a) and c)). The term 

1 r-' 

« = e-^^^ , 5, = - / \q{p,E)\dp (11) 

accomits for tunneling through the barrier, 

S4, = avgT{^+iS,)-S,\og\S,\+S, (12) 

is a phase correction, and Sg — below the barrier. Deep 
below the barrier, tunneling can be neglected and the 
simple semiclassical single well quantization is recovered 
(sec also @). 

Above the barrier, the inner turning points p^'' , p^'' 
turn into a complex conjugate pair and different contin- 
uations of the semiclassical quantization have been sug- 
gested [Hill [11]. Following [13 we replace these turn- 
ing points by the momentum at the barrier Pbarr in the 
formulas for Si^r, modify the tunneling integral Se as 

^ kp^:^ ~ pf) - - q{p,E)dp (13) 

and introduce a non- vanishing action integral 

/•Pbarr /"Pbarr 

Se = / q{p, E)dp+ / q{p, E) dp. (14) 

The combined semiclassical approximation is continuous 
if the energy varies across the barrier (from region II to III 
in Fig. [5]) and continuously approaches the simple version 
with only two turning points p j and in region III 
high above the barrier. 

Figure [3] shows the semiclassical many particle energy 
eigenvalues in dependence on the parameter e in the su- 
percritical region for = 10 particles {g = —3/Ns, 
V = 1, h — 1). Also here one observes an almost pre- 
cise agreement with the exact eigenvalues and the net of 
avoided level crossings in all details. In particular the 
level distances at the avoided crossings are reproduced 
and allow furthermore a direct semiclassical evaluation. 
Figure [S] shows both exact and semiclassical eigenvalues 
in dependence on e for subcritical interaction for only 
N = 2 particles. Even for that small particle number 
the semiclassical eigenvalues approximate the exact ones 
quite well. The deviation between the semiclassical and 
exact quantum eigenvalues decreases with increasing par- 
ticle number A^. 

For a more quantitative comparison, the exact and 
semiclassical eigenvalues are listed in Table U for A^ = 20 
particles and selected e-values. Here the relative error is 
only about 5 • 10~^. 
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FIG. 6: (Color online) Exact (- -), En and semiclassical ( ) 
many particle energies, E^'^, as a function of the onsite energy 
e in the subcritical region {g = — 0.9/A^s) for v = 1 and N = 2 
particles {h= 1) . 
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TABLE 1: Exact quantum -En=o...20 and semiclassical eigen- 
values for £ = 0, 0.5, 1.0, 1.5 for = 1 and A = 20 
particles (h — 1) in the supercritical region {g — ~3/Ns). 



C. Eigenfunctions 

The mean-field approximation allows also a semiclas- 
sical construction of the eigenstates = £'„|^„) 
of the Bosc-Hubbard Haniiltonian ^ resp. In the 
quantum case, the main interest may concentrate on 
the population imbalance of these states, i.e. the p- 
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representation 



p=-N:2:N 



(15) 



where p runs from — to in steps of 2. Based on the 
(classical) mean-field dynamics, we have to construct a 
semiclassical approximation in momentum space, which 
is discussed in some detail in [s^l- The purely classical 
momentum probability distribution is easily calculated 
as w^^ip) {2T\dn/dq\)-\ where T is the period of 
oscillation. Note that the factor 2 arises from the two 
classical contributions, i.e. the direct path and the path 
once reflected at the opposite turning point. For our 
mean-field Hamiltonian ^ this leads to 



w^^ip) = C v^{K 



(16) 

in the classically allowed region, where C = 1/2T takes 
care of the normalization. The so-called primitive semi- 
classical wavefunction includes interference of the two 
classical paths: 

\K'{p)\'' ^2\wM\ cos^ {Sip)/h~7:/A) (17) 

where S{p) = S{p,E^) is the classical action for an 
energy equal to the semiclassical eigenenergy of 
state number n, i.e. the oriented momentum-integral over 
q{p)=q{p,E^^-) 



Sip) 



q{p')dp' 



[P-P-] 



if p_ lies on the lower potential curve U- or 
S{p)^ r q{p')dp'. 



(18) 



(19) 




FIG. 7: (Color online) Momentum "potentials" U±{p) and 
- plotted at the energy En - exact (blue circles), primitive 
semiclassical (green) and uniform semiclassical (red crosses) 
wavefunctions |^n(p)|^ of the lowest eigenstate n = for = 
14 particles [g = -0.6 /Ns, e = 0.6, u = 1 and = 1). The 
solid lines are drawn to guide the eye. 




FIG. 8: (Color online) Momentum "potentials" U±{p) and 
- plotted at the energy E„ - exact (blue circles), primitive 
semiclassical (green) and uniform semiclassical (red crosses) 
wavefunctions |^'„(p)p of eigenstate n = 2 for A'^ = 14 parti- 
cles {g = —0.9/Ns, e = 0, V = 1 and h = 1). The solid lines 
are drawn to guide the eye. 



if p_ lies on U-^-. In the classical forbidden region q{p) is 
complex valued and we can use the usual complex con- 
tinuation 1291 



\K'{p)\' 



1 



<'(p)cxp(-2i5(p)/ft)| 



where 



S{P) 



T Jp q{p) dp, p < p_, j5_ on U± 
T /p+ 9(p) dp, P> P+, p+ on U± 



(20) 



(21) 



Note that these distributions should be renormalized to 
unity. 

The divergence at the classical turning points p± is 
finally cured by a uniform semiclassical approximation 
(see e.g. [l^). Here the different turning point scenarios 
discussed above must be distinguished. In the following we 
only state the result if p± both lie on the lower potential 



curve U^{p). A map ping onto harmonic oscillator wave 
functions then yields [2^ 



u';^ {p) ^2n + l- e H„ (^^ exp (-^ 



(22) 

where H„ denotes the Hermite polynomial of ?ith order 
and ^ is determined by 



1 



1 



with ^0 ~ \/2n + 1. 

Up to now, the semiclassical momentum variable p 
could be treated as continuous in the mean-field approx- 
imation, whereas in the quantum system, p is a dis- 
crete variable, p = —N, —N + 2, . . . , +N. Semiclassi- 
cally, this is a consequence of the even symmetry and the 
TT-pcriodicity of the mean-field Hamiltonian ^ in the co- 
ordinate q. As in Fourier transformation, this allows only 
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even or odd integer values of p. 

The final uniform semiclassical wave functions in mo- 
mentum space are therefore given by (|22|) at p = 
-N,...,N, normalized as Ep=-jV:2:Ar l*n b)P = 1- 

Figures [7] and [8] show a comparison of the primitive 
semiclassical approximation (normalized to fit the cen- 
tral maximum) and the uniform one with exact quantum 
results, both in the subcritical region for = 14 parti- 
cles. Shown is the ground state n = for a biased Bose- 
Hubbard model (e = 0.6) and the third excited state 
n = 2 for a symmetric one (e = 0). As expected, the 
quantum distributions mainly populate the classically al- 
lowed region inside the "potential" curves U±{p) and are 
very well approximated by the primitive semiclassical dis- 
tributions. In particular, the imiform approximation is 
almost indistinguishable from the exact values. 

IV. CONCLUSION 

It is demonstrated for a two-mode Bose-Hubbard 
model, that the mean-field approximation can be used to 
reconstruct approximately the individual eigenvalues in 
a semiclassical Bohr-Sommerfeld (or EBK) manner with 
astounding accuracy even for a small number of particles. 
The same holds for the primitive semiclassical approxi- 
mation of corresponding eigenstates which was shown for 
the subcritical case. Furthermore the possibility of a uni- 
form approximation was demonstrated for a special case. 
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